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By means of parallel tempering Monte Carlo simulations we find strong evidence for a finite- 
temperature spin-glass transition in a system of diluted classical Heisenberg dipoles randomly placed 
on the sites of a simple cubic lattice. We perform a finite-size scaling analysis of the spin-glass 
susceptibility, Xsg, and spin-glass correlation length, £l. For the available system sizes that can be 
successfully equilibrated, the crossing of £l/L versus temperature is strongly affected by corrections 
to scaling and possibly by a short-length scale ferromagnetic spin blocking. Similarly to many 
studies of different three dimensional spin-glass systems, we do not find a crossing of the spin glass 
order parameter Binder ratio. 



INTRODUCTION 



Our most formal current theoretical understanding of 
the spin-glass (SG) phase is based on the replica symme- 
try breaking (RSB) picture set by the Parisi solution 1 ' 2 of 
the infinite-dimensional Sherrington-Kirkpatrick model. '' 
As the upper critical dimension (UCD) of SG models is 
large (g?ucd=6),' such mean-field description is likely to 
be unsuitable to understand the physics of real materials 
exhibiting glassy behavior. An alternative description of 
the SG phase in finite dimension is given by the phe- 
nomenological droplet picture, ' which has been found 
to possibly characterize better three dimensional (3D) 
SG. 6 However, it remains an open debate what is the 
proper theory describing SG systems in real (finite) di- 
mensions. Most of our knowledge about the properties of 
3D SG models has been assembled thanks to years of ex- 
tensive numerical simulations. Unfortunatelly, the slow 
relaxation characterizing spin glass systems makes the 
numerical studies very difficult. Furthermore, the largest 
fraction of the numerical work has so far concentrated 
on the Edwards- Anderson (EA) model 7 of n-component 
spins interacting via nearest-neighbor random exchange 
interaction, Jy, where both ferromagnetic or antiferro- 
magnetic couplings are present. The cases n=l and n=3 
refer to Ising and Heisenberg SG, respectively. The prob- 
ability distribution of the random bonds, P(Jij), is usu- 
ally taken to be Gaussian or bimodal. 4 ' 8 

A great part of the numerical studies of SG models 
has been devoted to the minimal EA model, the one- 
component Ising SG. Due to severe technical difficulties, 
only very limited range of system sizes was accessible 
in the early simulations, while scaling corrections in SG 
systems are large. As a result, the existence of a fi- 
nite temperature SG transition in 3D Ising SG model 
had remained under debate for a long time. 9 ' 10 ' 11 ' 12 ' 13 ' 14 
The early MC studies strongly supported the finite- 
temperature SG transition, but zero-temperature tran- 



sition could not be definitely ruled out. 9 ' 10 ' 15 Only quite 
recently, in the course of large-scale Monte Carlo studies, 
has the existence of a thermodynamic phase transition in 
the Ising case been seemingly firmly established 11 and, 
perhaps most satisfactorily, universality among systems 
with different bond distributions been confirmed. 12 ' 13 

The case of the Heisenberg SG still remains somewhat 
more controversial than of the Ising SG. Originally, it 
was believed that the lower critical dimension (LCD) for 
the Heisenberg SG is c?lcd > 3, 16 and that the small 
anisotropies present in the real system are responsible 
for the SG behavior observed in experiments. 10 ' 1 ' While 
compelling, this suggestion has some difficulties since no 
crossover from Heisenberg to Ising SG universality class 
caused by a weak anisotropy has ever been observed in 
experiments. 18 ' 19 It was suggested that in the Heisen- 
berg EA SG, a finite-temperature transition occurs in 
the chiral sector 18 ' 19 while a SG transition in the spin 
sector occurs at zero temperature. 19 ' 20 ' 21 The chirality 
is a multi-spin variable representing the handedness of 
the noncolinear or noncoplanar spin structures. 18 ' 19 It 
has been proposed that the SG phase in Heisenberg SG 
materials is caused by a spin-chirality coupling induced 
by small anisotropies. 19 Later simulations indicated the 
existence of a nonzero-temperature SG transition. 22 ' 23 
The most recent work on the 3D Heisenberg SG sug- 
gests that the SG transition may be decoupled from the 
chiral glass (CG) transition, occurring at slightly higher 
temperature, ' 2o or that a common transition tempera- 
ture may exist. 21 ' In two dimensions, recent defect-wall 
renormalization group calculations suggest that both the 
SG and CG transitions occur at zero temperature, but 
that the two transitions are decoupled. 2 '' 28 ' 29 

On the experimental side there exist some SG ma- 
terials with a strong Ising-like uniaxial anisotropy, but 
the majority of experimental SG studies focus on nearly 
isotropic, Heisenberg-like systems. A well studied Ising 
SG material is Feo.sMno.sTiOa, 30 ' 31 as other examples 
of an Ising SG more recently found, Euo.sBao.sMnOa 
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and Cuo.5Coo.5Cl2-FeCl3 ?, ' ! ' !4 can be mentioned. Con- 
sidering Feo.sMno.sTiOa, 30 ' 31 the leading coupling is a 
nearest-neighbor exchange, and the compounds, FeTi03 
and MnTiOs, are antiferromagnets. In both cases, 
the nearest-neighbor exchange interactions within the 
hexagonal layers is antiferromagnetic. The magnitude 
of the intralayer coupling is substantially larger than 
the interlayer coupling. The SG nature of the mix- 
ture, Feo.5Mno.5Ti03, originates from the fact that the 
coupling between the layers is ferromagnetic in FeTi03 
but antiferromagnetic in MnTiOs; 30 hence, in the mix- 
ture, random frustration occurs. For the Heisenberg 
case, some short-range SG compounds are also avail- 
able, for example insulating EuzSri-aS. 4 ' 35 In the Eu- 
rich case, this material is a ferromagnet; the nearest- 
neighbor exchange interaction between Eu 2+ ions is fer- 
romagnetic and the next-nearest-neighbor exchange is 
weaker and antiferromagnetic. 4 ' 36 When magnetic Eu is 
randomly substituted with nonmagnetic Sr, a random 
frustration of the ferromagnetic and antiferromagnetic 
bonds arises. But the most often studied SG are nearly 
isotropic (Heisenberg) metallic systems interacting via 
a long-range Ruderman-Kittel-Kasuya-Yoshida (RKKY) 
interaction between localized magnetic moments medi- 
ated by conduction electrons. In this category, the classi- 
cal systems are alloys of noble metals such as Ag, Au, Cu 
or Pt, doped with a transition metal, such as Fe or Mn, 
often labeled as canonical SGs. In the large r limit, where 
r is the distance separating the magnetic moments, the 
RKKY interaction varies with r as cos(2fcp7 , )/r 3 , where 
hp is the Fermi wavevector. 

An another class of SG materials consists of spa- 
tially disordered magnetic dipoles. The dipolar interac- 
tion has either ferromagnetic or antiferromagnetic char- 
acter depending on the relative position of the inter- 
acting dipoles. In the presence of positional disorder, 
this gives rise to random frustration and a SG phase 
at low temperature and sufficiently high level of disor- 
der is expected. 1 ' A number of dipolar Ising SG ma- 
terials have been identified and related models have 
been studied numerically. With the aim of modeling 
nanosized magnetic particles dispersed in a frozen non- 
magnetic solvent, 40 systems of Ising dipoles on fully 
occupied 11 and diluted 12 (with x=35% and x=50% oc- 
cupancy) simple cubic (SC) lattice with randomly ori- 
ented easy axes have been simulated. In three dimen- 
sions, a spin-glass transition has been identified, both in 
the diluted 12 and undiluted case. 11 A well known physi- 
cal realization of a diluted dipolar Ising model and dipo- 
lar SG is LiHoaYi-aTM. 43 ' ' ' 4e Early on, some authors 
suggested the existence of an exotic anti-glass phase at 
very low concentration, x, 43 ' 47 '* 18 or questioned the exis- 
tence of SG transition in LiHo :c Yi_ 2 ;F4 altogether. 19 ' 00 
As well, early numerical studies of diluted Ising dipoles 
on SC lattice 51 and for a lattice geometry corresponding 
to LiHoaYi-aT^''-' did not find a spin- glass transition in 
diluted dipolar Ising systems. 01 ' 52 A more recent work, 
however, reports a spin-glass phase in a model approx- 



imating LiHo^Yi-zF^ 46 As in previous work, 52 cross- 
ing of the spin-glass Binder ratio plots was not found in 
Ref. [46], but a finite-size scaling of the spin- glass corre- 
lation length provided a compelling evidence for a ther- 
modynamical phase transition. 40 Apparently, the correc- 
tions to scaling are large for the sizes of dipolar systems 
studied, being very pronounced in the Binder ratio while 
the SG correlation length is somewhat less affected. 4 ' 1 
Very recent MC simulations of a site-diluted SC lattice of 
Ising spins coupled via a long-range dipolar interactions 
also found a finite-temperature spin-glass transition, but 
with different value of the correlation length exponent 
v=Q. 95, 53 compared with v=1.3 reported in Ref. [46]. 

The case of dipolar Heisenberg SG is an obvious ex- 
tension of the SG phenomenology reviewed above. In 
the presence of spatial disorder, the off-diagonal terms in 
the dipolar interaction destroy the rotational symmetry 
of the ground state. Thus, the dipolar Heisenberg SG 
is expected to be in the Ising universality class. 17 ' 54 ' 55 In 
this context, it would be interesting to study a three com- 
ponent (n=3, Heisenberg) SG system where anisotropic 
long-range interactions, i.e. dipolar interactions, domi- 
nate. Finding in such system critical exponents that are 
consistent with the exponents of Ising SG universality 
class would further confirm universality in spin glasses 
and boost our confidence in our largely numerically-based 
understanding of real spin glass systems. 

Experimentally, a diluted dipolar SG can be realized by 
sufficiently diluting magnetic dipoles with a nonmagnetic 
substituent, to the point that a short-range exchange in- 
teraction becomes insignificant, and long-range dipolar 
interaction dominates. The best candidate materials are 
compounds containing rare earth magnetic ions, as due 
to the screening of the partially filled 4/ shell by outer 
shells, the exchange interaction is relatively weak among 
rare earths, while their magnetic moments can be large. 

It was mentioned above that Eu^Sri-^S, at concentra- 
tion x ~ 0.5, is an example of a short-range Heisenberg 
SG, where the SG freezing is driven by the frustrated 
nearest-neighbor and next-nearest-neighbor exchange in- 
teractions. The 4/ electrons of the rare earth Eu 2+ ion 
give an a S 7 / 2 ground state with a sizable magnetic mo- 
ment of 7/iB- 36 Below the percolation threshold, i.e. x c = 
0.136 for the face centered cubic lattice (FCC) with first- 
and second-nearest-neighbor interactions, 56 the existence 
of a dipolar SG in Eu^Sri-^S was suggested. 5. To explain 
the two maxima in the ac susceptibility, at temperatures 
of order of 100 mK and 10 mK, an interplay of dipolar 
freezing and blocking of small clusters was proposed. " 
But the authors of subsequent studies ' 8 suggested that 
the features in the ac susceptibility of Eu^Sri-^S at con- 
centration 0.05 < x < 0.13 should be interpreted as a 
spin blocking and not a dipolar SG freezing. Neverthe- 
less, it was also proposed that maybe at lower concen- 
tration, x, a dipolar SG freezing in this material could 
be studied. 58 It should be reminded that the existence 
of a dipolar SG at high dilution in LiHo a; Yi_ a ;F4 has 
also been much questioned over the years. 43 ' 44 ' 5 ' 49 ' 50 
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Notwithstanding that recent experimental work reports 
a SG transition in LiHo a; Yi_ I :F4(x=0.045), 15 it is in- 
teresting to ask if difficulties in establishing the exis- 
tence of a dipolar SG in strongly diluted Eu^Sri-^S 57,58 
and in LiHo^Yi-^F^ i.e. the suggested anti-glass phase 
at a;=0.045 43 > 47 ' 48 or absence of signatures of SG tran- 
sition both at £=0.045 and x=0.165, 19, °° are related. 
A clearer picture of the formation of relatively large 
spin blocks below the percolation threshold that collec- 
tively form a cluster SG phase was obtained from exper- 
iments on Eu^Cai-aBg.' 9 Clusters with magnetic mo- 
ments [i ~ 260^tB were observed and the transition tem- 
perature separating the cluster glass and paramagnetic 
phases, for Eu concentrations between around x=0.1 and 
x=0.3, was measured to be of order of 2 K.° 9 

Promising candidates for diluted dipolar Heisenberg 
SGs can be found among gadolinium compounds. Gd 3+ 
ion has a half-filled 4/-shell. The ground state mani- 
fold is 8 5*7/2 and with little orbital momentum contribu- 
tion. Gd 3+ is therfore a good approximation of a classi- 
cal Heisenberg spin. Good example of materials that can 
be considered as candidates for diluted dipolar Heisen- 
berg SGs are (Gd x Yi_ K ) 2 Ti 2 7 and (Gd^Yi-^S^Or. 
Gd2Ti207 and Gd2Sn2C>7 are strongly frustrated Heisen- 
berg pyrochlore antiferromagnets. 60 ' 61 ' 62 While the 
Curie- Weiss temperature is about Oqw ~ —10 K, due to 
frustration, both compounds remain disordered down to 
T =1 K. 60,63 Theoretically, the extensive ground-state 
degeneracy in the pyrochlore nearest-neighbor Heisen- 
berg antiferromagnet prevents ordering down to zero 
temperature, 64,65 and the low temperature order in the 
aforementioned materials is induced by other, weaker in- 
teractions that are specific to each of these compounds. 
One of the interaction at play below 1 K is the dipo- 
lar interaction. Indeed, in the case of Gd2Sn2C>7, the 
spin configuration in the ordered state was found 66 to be 
the ground state of the pyrochlore antiferromagnet with 
dipolar interactions.'" In the case of Gd2Ti 2 07, other 
interactions, like further nearest-neighbor exchange, are 
likely at play. Indeed, below the first phase transition at 
1 K, there is another one at 0.7 K, 63 with both phases 
ordered with propagation vector k = [535] . 68,69 

Another candidate material, also well studied in the 
context of geometrical frustration, is the GdsGasO^ 
garnet (GGG). The cubic lattice structure of this frus- 
trated Heisenberg antiferromagnet consist of two inter- 
penetrating sublattices of corner-sharing triangles. While 
the Curie- Weiss temperature is 6*cw ~-2 K, the frus- 
tration postpones ordering down to 0.18 K. 70 The rich 
low temperature physics of GGG is still not fully under- 
stood, but some insight has been recently gained from 
dynamic magnetization studies, revealing that in the 
low temperature phase there is long-range order coex- 
isting with both spin liquid 11,1 1 and spin glass 73 behav- 
ior. In the framework of Gaussian mean-field theory, it 
was shown that the dipolar interaction plays an impor- 
tant role in the ordering in GGG, 74 ' 75 and that the neu- 
tron scattering data'" can be reproduced with a proper 



treatment of the dipolar interactions.' 4 '' 5 Analogously 
to (GdxYi_x)2Ti 2 7 and (Gd K Yi-x)2Sn 2 7j at suffi- 
cient dilution, (GdaYi-^GasO^ may be expected to 
exhibit, at low temperature, a dipolar SG phase. 

In the studies motivated by experiments indicating a 
SG freezing in frozen ferrofluids, 40 some evidence for a 
SG freezing in a system of dense amorphous Heisenberg 
and XY spins coupled by long-range dipolar interactions 
was obtained from molecular dynamics simulations. 38,39 
But no systematic investigation of the thermodynamic 
nature of the freezing was at that time really possible. 

In the present work, in anticipation of eventual experi- 
mental studies of dipolar SG, e.g. diluted Gd compounds, 
or further work on Eu^Sri-^S or Eu^Cai-^Be, we per- 
form numerical studies of the SG transition in a diluted 
dipolar Heisenberg model. At high dilution, the lattice 
structure should be irrelevant and data obtained for dif- 
ferent systems should be comparable. 76 Here we consider 
the simplest possible geometry where we study dipoles 
randomly placed at the sites of SC cubic lattice. We pro- 
vide Monte Carlo data that supports the scenario that, 
at low dipole concentration, the diluted dipolar Heisen- 
berg model displays an equilibrium phase transition to a 
SG phase. We calculate the critical exponents v and 77 
for the SG transition in the model studied. The derived 
exponents do not match experimental or Monte Carlo ex- 
ponents, neither for Heisenberg nor Ising SG. This may 
be because of important scaling corrections and severe 
restriction of the system sizes that we were able to study 
due to the computationally expensive summation of long- 
range dipole-dipole interaction and very slow equilibra- 
tion. 

The rest of the paper is organized as follows. In Sec- 
tion II, we define the model and MC method employed. 
In Section III, we introduce the observables calculated in 
the simulation. In Section IV, we present and discuss our 
results. Our conclusions are in Section V. Some techni- 
cal details are discussed in Appendices. In Appendix A, 
we examine the temperature and system size dependence 
of the magnetization and staggered magnetization of the 
model studied. In Appendix B, we discuss the issue of 
the self-interaction term that must be taken into account 
when periodic boundary conditions are imposed. In Ap- 
pendix C, we discuss the Ewald summation technique. In 
Appendices D and E, we discuss the overrelaxation and 
heatbath algorithms, respectively. 



II. MODEL AND METHOD 
A. Model 

We consider a system that consists of classical three- 
component (n=3, Heisenberg) dipoles that are free to 
point in any direction. The dipoles are randomly dis- 
tributed on the sites of a 3D simple cubic (SC) lattice. 
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The Hamiltonian is of the form 

,2 _3_M 



1 



i J 



r- 



S»( ri )S v ( rj ), (1) 



where S^fa) (fi — x, y, z) denote Cartesian components 
of classical spin vectors, S(r,), which are of unit length, 
|<S(Tj)| = 1. The energy scale of dipolar interactions is set 

2 

by ed = , where /i is the magnetic moment of the spin 
S(ri), a is the lattice constant and /xo denotes vacuum 
permeability. Below, ed/k-Q, where fcs is the Boltzmann 
constant, is conveniently used as a unit of temperature. 
The summation is carried over all occupied lattice sites 
and over the vector components of the spin, (1,1/ = x, 
y and z. The factor 1/2 is included to correct for dou- 
ble counting of dipole pairs, and rj are the positions 
of ions labeled i and j, respectively, and their distance, 
\rij\ = \rj — n\, is measured in units of nearest-neighbor 
distance, a. We use periodic boundary conditions. In the 
case of long-range interactions, this means that to calcu- 
late a pairwise interaction, we must sum over an infinite 
array or dipole images replicated with a periodicity set 
by the dimensions of the simulation box. Therefore, it is 
convenient to consider the interaction constant for spins 
i and j as a 3 by 3 matrix, Ly. The matrix elements 
of Lij are denoted , 



and for dipoles separated by a 



vector Vij , are given by the sum, 

m = ^ \m + "i 2 - 3 (m +™f fa + »r m (2) 



Vectors n are of the form n = L (kx + ly + mz), where 
k, I, m are integers and x, y and i are the unit vectors 
pointing in the directions of primitive translation vec- 
tors of the SC lattice. L is an integer expressing the 
size of the simulation cell in units of the lattice constant, 
a. Note that, in a simulation, care must be taken to 
correctly include the so-called self- interaction term, La, 
(see Appendix B). The self-interaction term originates 
from interaction of spin with its periodic images repli- 
cated outside the simulation cell. The lattice summation 
(2) is performed using the Ewald technique;' 7 ' 78 ' 79 ' 80 the 
details of which are given in Appendix C. The calculated 
Ewald sums correspond to a summation over a long cylin- 
der, such that the demagnetization field is zero. Using 
interaction constants defined in Eq. (2), the Hamiltonian 
can then be written in the form 
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H = -e d J2S(r i )L ij S(r j ). 



(3) 



The summation in Hamiltonian (3) includes only the 
spins enclosed in the simulation cell, while the presence of 
spins outside the simulation box is approximated by pe- 
riodic images of the spins within the simulation cell, and 
this effect is included in the interaction constants, given 
by the matrix calculated via the Ewald method. Note 
that, unlike in Eq. (1), the self-interaction terms, % = j, 
are included in the summation (3). 



B. Method 

Our simulations employ the standard single spin-flip 
Metropolis Monte Carlo algorithm with parallel temper- 
ing (PT). 81 ' 82 PT was found to significantly speed up 
equilibration in slowly relaxing systems. 81 ' 82 In this tech- 
nique, one simultaneously simulates a number, Nt, of 
thermal replicas - copies of the system with the same 
spatial disorder, but at different temperatures. In each 
thermal replica, the simulation begins from a different 
random initial spin configuration. At every 10 local up- 
date sweeps, each consisting of N single spin updates, 
where N is the number of spins in the system, a config- 
uration swap among thermal replicas is attempted with 
the acceptance probability preserving the detailed bal- 
ance condition. The frequency of tempering is chosen to 
balance the following two factors. As thermal temper- 
ing is computationally inexpensive, it is desirable to per- 
form replica swap attempts often, to promote traveling 
of the replicas along the temperature axis. But, on the 
other hand, after a parallel tempering induced configura- 
tion exchange, a sufficient number of local update sweeps 
has to be performed to let the new configuration evolve 
at the given temperature. If a subsequent tempering 
is attempted too soon, the two configurations could be 
swapped back, and in that case no progress in the relax- 
ation of a state trapped in a local energy minimum would 
have been made. The number of thermal replicas, Nt, 
and simulated temperatures, T a , where a = 1, . . . , Nt, 
are chosen to yield a sufficiently high and temperature 
independent PT configuration swap acceptance rate, i.e. 
not less than 50%. A uniform with respect to tempera- 
ture, T a , PT acceptance rate is achieved by choosing T a 
to satisfy the formula 23 



(T Q - T a _ x )lT a = 1/^/CvN, 



(4) 



were N denotes the number of dipoles. The specific heat 
per spin, Cy, used in Eq. (4) was measured in prelimi- 
nary simulations of the smallest system sizes, with uni- 
formly distributed temperatures. 

The Metropolis single-spin moves are attempted within 
a temperature-dependent solid angle, where the angle is 
self-consistently chosen such that the acceptance rate of 
MC single spin move is close to 50%. To carry out a spin 
move, we choose a coordinate system with the z axis 
along the current spin direction, and randomly choose a 
polar angle, 9, and azimuthal angle, 4>. In order to obtain 
a uniform distribution of random points on a unit sphere 
one needs to draw (ft and z — cos(6>) from a uniform prob- 
ability distribution, such that (ft £ (0, 2tt) and z S (—1,1). 
Here, to maintain the desired acceptance rate, the move 
is restricted to a limiting angle, 9 ma , x , relative to the ini- 
tial spin direction; hence, the choice of z is restricted 
to z e (1 - z max , 1), where z max = cos(6> max ). To ob- 
tain z m ax such that the acceptance rate, p acc , is 50%, 
during each 100 MCS p acc is measured, and afterwards 
z mal is adjusted. If p acc is lower than 0.5, z max should 
be decreased; in the opposite case, when p acc > 0.5, 
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L 


iVdip 


iVsamp 


iVcq 


A'prod 


N T 


7m in 




z=0.0625 


8 


32 


5000 


5-10 5 


5-10 5 


16 


0.05 


0.1763 


10 


62,63 


2000, 2000 


2-10 6 


10 6 


16 


0.05 


0.1763 


12 


108 


1200 


10 7 


10 6 


16 


0.05 


0.1763 


14 


172 


1000 


10 7 


10 6 


16 


0.05 


0.1763 


z=0.125 


6 


27 


3000 


5-10 5 


5-10 5 


16 


0.0750 


0.2869 


8 


64 


2000 


2-10 6 


10 6 


16 


0.0750 


0.2869 


10 


125 


1500 


2-10 6 


10 6 


16 


0.0800 


0.2811 


12 


216 


1000(+200) 


5-10 6 (2-10 7 ) 


10 6 


16 


0.0850 


0.2787 



Table I: Parameters of the Monte Carlo simulations for two 
dipole concentrations, x. L is the linear size of the simulation 
box; iVdip is the number of denotes number 

of disorder samples. N eq and A'prod are the number of MCS in 
the equilibration and measurement phase of the simulation, 
respectively. Nt is the number of thermal replica and T m i n , 
Tmax are the lowest and highest temperatures in PT scheme. 
For L=10, a;=0.0625, to obtain the desired x, two numbers of 
dipoles were simulated, and the disorder average was taken 
over the results for both A'dip = 62 and A'dip = 63. For L=12, 
1=0.125 the numbers given in round brackets pertain to a 
subset of disorder replicas simulated longer, to monitor equi- 
libration; the long equilibration time results for these replicas 
were included in the disorder averaging. 

z max should be increased, while for p acc = 0.5, z max 
does not change. Such update of z max can be obtained 
when multiplying the current value of z m ax, , by 

2pacc; hence, a new value of z max is calculated according 
to the formula zl™' = 2p acc Zmax\ with the restriction 
Zm™' € (0.001, 2). After choosing <p and 9, that is a new 
spin direction in the coordinates relative to the initial 
spin direction, a transformation to the global coordinate 
system is performed. 

We simulated two dipole concentrations, x = 1/16 = 
0.0625 and x = 1/8 = 0.125, and for each concentration 
we considered 4 system sizes varying between around 30 
and 200 dipoles, which is the largest size that we were 
able to equilibrate. To perform the necessary disorder 
average (see Section III), we considered at least 1000 dis- 
order samples. The parameters of the simulations are 
collected in Table I. To generate results reported here, 
we used in total around 3 • 10 5 hours (^35 years) of CPU 
time on AMD Opteron, 2.6 GHz. The statistical error is 
based on disorder sampling fluctuation and is calculated 
using the standard jackknife method. 83 ' 8 ' 85 

To reduce the number of performed lattice sums, for 
each lattice site k we calculate the local interaction field, 

H k = ^2L kj S j} (5) 

and update it only when the spin change is accepted. 
Having H k available, the computational complexity of 
calculating the energy change when a single spin is 



moved, which is needed to test if the spin move is ac- 
cepted, is of order of a small constant number of arith- 
metic operations, 0(1). As the field, Hk, is updated only 
if the spin move is accepted, the computational cost of 
updating the local field, Hk, which is 0{N), is avoided 
if the spin move is rejected. In consequence the compu- 
tational cost of rejected Metropolis spin updates is neg- 
ligible. 



In was reported that the autocorrelation time can 
be substantially decreased in simulations of Heisenberg 
SG by performing computationally inexpensive overre- 
laxation (microcanonical) spin updates. 86 ' 87 Overrelax- 
ation updates are zero-energy spin moves that consist of 
180 degrees rotation of the spin around the local molecu- 
lar field. Including overrelaxation moves in simulations of 
short-range Heisenberg spin systems is beneficial because 
overrelaxation moves are faster than Metropolis updates. 
Here, in the case of long-range interaction, the compu- 
tational cost of overrelaxation moves would not be less 
then the cost of Metropolis spin flips, as most of the time 
is spent on updating the local interaction field (5), and 
the local interaction field has to be updated both after an 
accepted Metropolis spin flip and after an overrelaxation 
move. Furthermore, it is worth to note that in the gen- 
eral case of non-cubic geometry with periodic boundary 
conditions, where the self- interaction term is present (see 
Appendix B), an overrelaxation move does not preserve 
the energy (see Appendix D). 



In the case of the nearest-neighbor Heisenberg SG 
model, it is more efficient to use the heatbath 
algorithm 88 ' 89 ' 90 for local spin updates. In the heatbath 
algorithm spins are individually connected to a heatbath, 
i.e. a new spin direction, which is independent from the 
previous direction, is drawn from the Boltzmann proba- 
bility distribution for a spin in the local molecular field, 
Hk- Hence, in contrast to the conventional Metropo- 
lis method, in the heatbath algorithm, the computa- 
tional cost of rejected spin moves is avoided. For the 
current problem, the benefit of using the heatbath algo- 
rithm would not be high because the computational cost 
of rejected spin update attempts is negligible in com- 
parison with the computational cost of accepted updates 
which require recalculating the lattice sums in Eq. (5). 
Also, similarly to the case of overrelaxation moves, if the 
geometry of the simulation cell is not cubic, the self- 
interaction term in the Hamiltonian introduced by the 
periodic boundary conditions makes the heatbath algo- 
rithm impractical (see Appendices B and E). So, while 
we are considering here a cubic system, in order to keep 
our method general and to obtain results that are the eas- 
iest to compare with possible future simulations with dif- 
ferent lattice geometries, we decided to not use the heat- 
bath algorithm nor overrelaxation moves in this work. 



III. PHYSICAL QUANTITIES 

In spin-glass systems, the order parameter can be de- 
fined as an overlap between two independent, identical 
copies of the system. In the case of 3D Heisenberg spins, 
the overlap can be calculated for 9 combinations of the 
vector components. We write 



(6) 



where fj,, v = x,y,z and where a and (3 denote differ- 
ent copies of the system with the same random disor- 
der and that are simulated simultaneously, but indepen- 
dently. The wave-vector-dependent SG order parameter 



is 



q(k)= /$><-(fc) 



(7) 



In the case of EA Heisenberg SG, in addi- 



tion to spin, a chirality variable 



18,19 



is also 



considered, 19 ' 20 ' 21,22,23 ' 24 ' 25 ' 26 and the chirality overlap 
order parameter and further quantities defined using 
this order parameter are calculated. Here, for a diluted 
dipolar Heisenberg SG, we do not consider chirality. The 
chirality cannot be easily defined in a diluted system. 
Moreover, note that, as discussed in the Introduction, 
the anisotropy of the dipolar interaction, in the presence 
of spatial disorder, brakes the rotational symmetry. 
Because of that, the chirality, if it was defined, cannot 
be decoupled from the spin. 

Traditionally, the finite-size scaling (FSS) analysis of 
SG simulation data has been based on the calculation of 
Binder ratios, 91 ' 92 ' 9 '' 1 which, for an n=3 Heisenberg SG, 
is defined as: 24 ' 26 




9 



[(g(Q) 4 )] 

[<9(0) 2 )] 2 



(8) 



where (. . .) denotes thermal averaging and [. . .] is a disor- 
der average. The numerical factors in Eq. (8) are chosen 
such that at T = oo, assuming Gaussian distribution of 
g(0), Ul = 0, and at T = 0, where q(0) is not fluctuating, 
Ul is 1. Being a dimensionless quantity, Ul is expected 
to display FSS properties described by 12 



U L = X{L l l"{T - T g )), 



(9) 



where the scaling function 1 ' 1 X is an analytic function of 
its argument, and v is the universal correlation length 
exponent, such that there is no system size dependence 
outside the argument of the scaling function. Many re- 
cent works report that in the case of disordered spin glass 
systems, a better FSS analysis can be achieved when con- 
sidering the finite-size SG correlation length, ^.i 1 . 23 . 46 . 96 
In the context of Ising SGs, it was suggested that U l may 
not cross due to a lack of unique ground state 96 or be- 
cause it is too noisy (see footnote Ref. [97]), 11 being a 



quantity that requires evaluation of a four-point corre- 
lation function, as opposite to £l, that is defined using 
a two-point correlation function. It was also observed 
that scaling corrections are larger for U l than for £l ■ ' ' 
It is likely that in the case of Heisenberg SG large scal- 
ing corrections are the leading factor behind the lack of 
crossing of the Binder ratios. To proceed, we define the 
SG susceptibility 8 ' 11 ' 2 ' 5 as 



XSG(k) = N [(q(k) 2 )] 



(10) 



Assuming an Ornstein-Zernike form for the SG 
susceptibility, 98 



XSG (fc)(xi/(|fcr + r 2 ), 



(11) 



where |fe| <C l/£. We define a finite-size SG correlation 
length, 11 ' 99 £l, via 



Xsg(O) 
2sin(fc mi „/2) \ 

V XSG (k m in) 



- 1 



1/2 



(12) 



The correlation length divided by the system dimen- 
sion, £,l/L, similarly to the Binder ratio, is a dimension- 
less quantity that is expected to scale according to the 
relation 11 ' 12 ' 23 ' 26 



U/L = Y{L X /»{T-T g )), 



(13) 



where Y is once again a scaling function. Hence, at a 
putative SG transition temperature, T g , £l/L is expected 
to be size independent. 



IV. MONTE CARLO RESULTS 

A system of Heisenberg dipoles on a fully occupied SC 
lattice orders antiferromagnetically. 100 ' 101 To rule out a 
long-range order in the simulated diluted systems, we 
calculated the magnetization, M, and the staggered mag- 
netization, Mgtag. Both M and M sta g are small and de- 
crease with increasing system size. This indicates that 
their nonzero value is a finite-size effect and not a result 
of long-range ordering. More detailed discussion of M 
and Mgtag is given in Appendix A. 

We plot in Fig. 1 the temperature dependence of the 
Binder ratio, Ul, for £=0.125 and 2=0.0625, for differ- 
ent system sizes. The Binder ratio curves do not cross; 
hence, they do not provide indication of a phase transi- 
tion. Also, in some studies of other models, a crossing of 
the Binder ratios was not found, while the scaling invari- 
ance of the finite-size correlation length was established, 
indicating a transition to a SG phase. The magnitude 
of scaling corrections is different for different observables 
and they are likely to be larger for Binder ratio than 
for correlation length. In the simulation of the Ising EA 
SG 11 ' 12 Ul does cross, but the scaling corrections are 
found to be larger for Ul than for ^/L. 11 In the case 
of the site diluted EA Ising SG, 102 where scaling correc- 
tions are large in comparison with other Ising SG models, 
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Figure 1: (color online). Binder ratios for x=0.0625 (left) and 
2=0.125 (right) as a function of temperature. 

£l /L plots are crossing with large shifts between system 
sizes, while the Ul curves do not cross, but merge at 
low temperature. A similar effect has been seen in the 
studies of diluted dipolar Ising SG' 16, ° 2 - Ul plots do not 
cross, but they have a tendency to merge at low T, while 
ih/L plots intersect. 4 ' 1 In the case of isotropic Heisenberg 
EA SG, 24,25,26 the behavior of spin and chirality Binder 
ratios differs, but neither show a crossing, while the cor- 
relation length shifts between system sizes shows that 
scaling corrections are large. It is worthwhile to note 
that the form of Binder ratio plots, characterized by a 
dip to a negative value in the proximity of T g , resembles 
the Binder ratio plots for chirality (and not spin) in the 
Heisenberg EA model, 24 ' 26 or the Binder ratio plots for 
spin in Heisenberg SG models in the presence of random 
anisotropy in three 1 " 1 and four 104 dimensions. 
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Figure 2: (color online). SG correlation length as a function 
of temperature, £=0.0625. 

Having discussed the temperature dependence of the 
Binder ratios, we now turn to the behavior of the SG 
correlation length, £,l{T). We show the plots of £,l/L vs 
T for various system sizes in Figs. 2 and 3. The curves do 
cross; but, for both concentrations, there are large shifts 
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Figure 3: (color online). SG correlation length as a function 
of temperature, £=0.125. 



between the intersection points for different system sizes. 
The large shifts between the intersection points were also 
found in the studies of the EA Heisenberg SG. 24,25,26 In 
the case of EA Heisenberg SG, a broad range of system 
sizes were studied and the shifts between the intersection 
points were systematically analyzed. 24,25,26 The limita- 
tion with our data, which are due to time consuming 
summation of long range interactions, prevent us from 
investigating large system sizes and, therefore, to per- 
form such an analysis. Because of a narrow range of 
available system sizes, the separation between the curves 
in the crossing region is small in comparison with the er- 
rorbars. Hence, the statistical uncertainty of locating the 
crossing points would be large. Indeed, looking at Figs. 2 
and 3 one realizes that by moving the curves within the 
error bars, the position of the crossing can be changed 
substantially. Also, the number of system sizes and, con- 
sequently, the number of intersection points is small. In 
our vs T data, the shifts between the intersection 

points for the smallest system sizes, consistently for both 
concentrations, are much smaller than the shift of the in- 
tersection point of the two largest system sizes studied. 
Such feature have not been found in simulations of EA 
Heisenberg SG. 24,25,26 A possible explanation for such 
behavior may be existence of short-range ferromagnetic 
correlations. Such a ferromagnetic spin blocking would 
especially strongly affect the data for the smallest system 
sizes, which possibly have the linear dimensions compa- 
rable with the length scale of the short-range ferromag- 
netic correlations. Spin blocking has been observed ex- 
perimentally in studies of diluted dipolar Heisenberg sys- 
tems EuajSri-zS 57,58 and Eu^Cai-^Be. '" Another argu- 
ment supporting short range ferromagnetic correlations 
scenario in the model studied herein is the finite-size mag- 
netization found for small system sizes (see Appendix A) . 

The scaling equation (13) is expected to be satisfied 
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Figure 4: (color online). Extended scaling of £l/L at 
x=0.0625 
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Figure 5: (color online). Extended scaling of £l/L at a;=0. 125 



only in close proximity of T g . To better describe the 
data at a larger distance from the critical point, Camp- 
bell et al. proposed a heuristic extended scaling scheme 
(ESS) 105 for £ L /L of the form: 



£l/L = Y[ (TL) 



1 - 



T 
JL 

T 



(14) 



and showed the improvement of the accuracy it provides 
in the case of a 2D Ising ferromagnet. Based on the as- 
sumption of a symmetric interaction distribution, P(Jij), 
they proposed, and tested numerically, an alternative 
scaling formula for the Ising EA spin glass, where t s/t in 
Eq. (14) is replaced with ( t s/t) 2 . Here we use the scaling 
formula of Eq. (14) because the bond distribution in the 
case of diluted dipoles is not symmetric. In a recent MC 
simulation of a diluted dipolar Ising SG, an ESS as given 
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Figure 6: (color online). Conventional scaling of £l/L with 
L 1/v (T - T g ); a~=0.0625 (left) and x=0.125 (right). 



by (14) was also found to describe the scaling of £,l/L 
better than if ( t s/t) 2 was used. 46 

We fit our data for £l/L to the scaling function (14) 
over the whole simulated temperature range, shown in 
Table I. The scaling function, Y, is approximated with a 
6th order polynomial, 



*•(*) = E 



(15) 



m=0 



where z = (TL) 
function, 



l/u 



(1 — t h/t). We define the penalty 



D= J2 (F(z)L/Z L -l)\ 



(16) 



MCdata 



that is minimized with respect to the parameters {a m }, 
T c and v. We obtain the values of the critical exponent 
v =1.16, v =1.09, and transition temperatures T g =0.074, 
T 9 =0.12 for x=0.0625 and £=0.125, respectively. The 
scaling collapse of the simulation data is shown in Fig. 4 
and Fig. 5. 

In Figs. 6 and 7, just for comparison, we present the 
results of the fitting to the conventional formula (13) and 
ESS with ( t s/t) . The fitting to the conventional formula 
(13), shown in Fig. 6, gives quite similar results to the 
ESS of Eq. (14). Apparently, the inaccuracy due to 
the small system sizes studied is larger here than the 
correction made by replacing Eq. (13) with Eq. (14). 

In the case of fitting to Eq. (14) with ( t b/t) , we obtain 
a visibly worse data collapse than when t s/t is used; the 
result of such a fit is shown in Fig. 7. 

We plot in Figs. 8 and 9 the SG susceptibility for 
fe = 0, xsg(0), of Eq. (10) for x=0.0625 and x=0.125, 
respectively. 

The SG susceptibility is expected to scale according to 
the ESS formula 105 



XSG = (TL) 2 "" Z ( (TL) 



T 



(17) 



We performed a fit following a procedure similar to the 
method used for the scaling fit of £l/L described in 
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Figure 7: (color online). Extended scaling of £l/L with 
(TL) 1/v (l-{Tg/ T f); a~=0.0625 (left) and z=0.125 (right). 
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Figure 8: (color online). SG susceptibility, :r=0.0625 



Eqs. (15) and (16). For x=0.0625 we obtain T g = 0.078, 
j/=1.25 and 77=1.45. For £=0.125 we get T g = 0.12, 
z/=1.18 and 77=1.35. The critical temperatures are con- 
sistent with those obtained from FSS of The values 
of the critical exponent v obtained here are slightly larger 
than v obtained from the scaling of £,l/L. The scaling 
collapse of xsg is plotted in Fig. 10 and 11 for x=0.0625 
and 2=0.125, respectively. 

In the dipolar Hamiltonian (1) off-diagonal terms, that 
couple different vector components of the dipolar mo- 
ment, are present. The off-diagonal terms destroy the 
rotational (0(3)) symmetry in an otherwise isotropic vec- 
tor spin system, and only a Z2 symmetry remains. It was 
suggested that such spatially disordered dipolar systems 
belong to the Ising universality class. 1 '' 51 Due to spatial 
disorder, the couplings, including the off-diagonal terms, 
are random, and the distribution of local freezing direc- 
tion in the SG phase remains uniform, unlike in a system 
with a global single-ion anisotropy (e. g. ~DSl term in 
the Hamiltonian) . With the uniform distribution of local 
freezing directions, a system is said to have a statistical 
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Figure 9: (color online). SG susceptibility, £=0.125 
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Figure 10: (color online). SG susceptibility scaling, £=0.0625 



rotational symmetry. 17 

The values of the critical exponents found in this work 
do not agree with either those from simulations of short- 
range Ising SG, v = 2.45, 77 = —0.375, 13 nor Heisenberg 
SG, v = 1.49, 77 = —0.19. 26 It is possible that our expo- 
nent v is consistent with v = 1.3 and v = 0.95 obtained 
for Ising diluted dipolar SG in Ref. [46] and Ref. [53] , re- 
spectively. Extracting SG critical exponents from simula- 
tions is difficult. Critical exponents for the Ising SG have 
been discussed for a long time 9 ' 10 ' 11 ' 12 ' 13 ' 14 , and proposed 
values were changing much with progress in development 
of simulation algorithms and computer hardware. Simi- 
larly to the early simulations of the Ising SG, 9 ' 10 our data 
very likely suffer from large scaling corrections. As our 
system sizes are small, one may want to compare our ex- 
ponent v with the results of simulations of the Ising SG 
performed for small system sizes, e.g. these in Ref. [9] 
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Figure 11: (color online). SG susceptibility scaling, £=0.125 




Figure 12: (color online). Snapshot of 200 equilibrated in- 
dependently spin configurations for L=12, at T=0.05 and 
x=0.0625. The alignment, i.e the scalar product, of the spins 
with the local freezing axes is indicated by the colors of ar- 
rows. 



{y = 1.2) or Ref. [10] {v = 1.3). There is a fair agree- 
ment in v but not in r\. The value of the exponent rj, from 
simulation 11 ' 12 ' 13 ' 26 and experiments 4 ' 31 ' 32 on many dif- 
ferent materials, for both Ising and Heisenberg SG is a 
small number, either positive or negative but not exceed- 
ing 0.5 in absolute value. Surprisingly, the value of ?y we 
obtain for the diluted dipolar Heisenberg SG, rj — 1.4, is 
much larger. 

Having discussed the question of universality class and 
commented on the expectation that, for diluted n=2> 
component dipoles, it should be Ising, it is interesting 
to ask whether such Ising structure is explicitly physi- 
cally manifest in the low temperature regime of the sys- 



tems studied above. We show in Fig. 12 a number of 
super imposed snapshots of the spin configurations for 
one disorder realization, in the low temperature phase, 
at T=0.05, and dipole concentration £=0.0625. The im- 
age contains the spin configurations of 200 replicas of 
the same disorder, each equilibrated independently, start- 
ing from different random initial configuration. The sys- 
tem size is L=12, which, at the concentration x=0.0625, 
gives 7V=108 ions. The parameters of the simulations, 
i.e. temperatures and number of equilibration sweeps, 
are given in Table I. In the case of isotropic Heisen- 
berg models, a low temperature phase has 0(3) rota- 
tional symmetry, and one expects the spin directions in 
replicas of the same disorder as explained above to be 
uniformly distributed. Here, due to the anisotropic char- 
acter of the dipolar interactions, a subset of the dipoles is 
characterized by a unique Ising local freezing direction. 
It is indicated by the fact that in the snapshots some 
dipoles have a strong tendency to point along a particu- 
lar local random direction, i.e the arrows can be enclosed 
by a circular conical surface with a small opening angle. 
Such inhomogeneous "random Ising structures" have also 
been observed in a model of diluted two-component 2D 
quadrupoles. For clarity, the alignment of spins with the 
local freezing directions, which is measured as an absolute 
value of the scalar product of a spin and the local freez- 
ing direction, is indicated by the color of the arrows. The 
local freezing direction vector is computed by summing 
all the spin vectors at a given site for the 200 disorder 
realizations in the following way. Starting from the sec- 
ond element in the sum, it is checked if adding another 
vector to the existing sum will increase or decrease the 
magnitude of the new sum. If adding the new element 
is to decrease the magnitude of the sum, the spin vector 
is added with a minus sign, such that the magnitude of 
the sum always increases. In this way we obtain a vec- 
tor that is pointing along the local freezing axis. Not all 
the sites are characterised by a local freezing direction. 
The arrows on the sites that do not have a local freezing 
direction create spherical structures. These dipoles have 
freedom to point in any direction in the low temperature 
phase. That means that these dipoles are strongly frus- 
trated and decoupled from the other dipoles. It is inter- 
esting to note that this behavior resembles the presence 
of "protected degrees of freedom" observed in gadolinium 
gallium garnet (GGG).' 2 The sites with a local freez- 
ing direction, in Fig. 12, seem to form small clusters. 
Possibility of ferromagnetic spin blocking was mentioned 
earlier as a potential explanation of a large shift of the 
correlation length crossing points for the largest system 
sizes relative to the crossing points for the smaller sizes. 
Formation of ferromagnetic spin blocks is also suggested 
by nonzero, but decreasing with system size, finite-size 
magnetization (see Appendix A). 

In the simulations of a SG system, it is of paramount 
importance to ensure equilibration of the system before 
the statistics for the measured observables is collected. 
As the quantity of foremost interest here is the correla- 



11 



0.18 

10" 



L=8 


T=0.05 


L=10 




— L=12 




L=14 


7^1 5 — 3 


x=0.0625 



10 10 

MCS 




10 10 
MCS 



Figure 13: (color online). Equilibration £=0.0625 (left) and 
2=0.125 (right). 



tion length, we assume that the system is equilibrated 
when the correlation length reaches a stationary state. 
We plot in Fig. 13 £,l/L vs the number of the equi- 
libration steps performed before the measurement was 
taken. The number of necessary equilibration steps in- 
creases very rapidly with the system size and, because 
of this, we were only able to equilibrate system sizes up 
to about 200 dipoles. We observe that the long-range 
dipolar Heisenberg SG takes longer time to equilibrate 
than the short-range 3D EA Heisenberg SG model. 24 ' 26 
A similar fact has been observed in the case of dipolar 
Ising SG. 16 



V. SUMMARY 

In conclusion, we have studied the spin-glass (SG) 
transition in a diluted dipolar Heisenberg model. From 
an analysis of the finite-size scaling of the SG correlation 
length, £l, we found an indication of a SG transition at 
a temperature T g =0.074, T s =0.12, and critical exponent 
v =1.16, v =1.09 for dipole concentrations x=0.0625 
and x=0.125, respectively. From finite-size scaling of the 
SG susceptibility, xsg> we obtained T g = 0.078, ^=1.25, 
77=1.45, and T g = 0.12, u=1.18, 77=1.35 for x=0.0625 and 
£=0.125, respectively. As in the isotropic Heisenberg SG, 
the Binder ratios, Ul, do not exhibit a crossing for dif- 
ferent system sizes. 

Our data support the scenario of ferromagnetic spin 
blocking. Short-range ferromagnetic correlations are in- 
dicated by a relatively large finite-size magnetization. 
Such short-range correlations would also explain unusual 
behavior of the SG correlation length The crossing 
points of the S,l/L vs T plots for the largest system sizes 
is shifted to much lower temperatures from the crossing 
points for the smaller system sizes. It may be caused by 
reaching a system size that is larger than the length scale 
of ferromagnetic clustering. Some indication of formation 
of frozen spin clusters can be also found from inspection 
of spin configuration snapshots. 

The long-range interactions, and hence the large num- 
ber of interacting spin pairs, give rise to a larger level of 
random frustration than in short-range (nearest neigh- 



bor) SG. Diluted dipolar SG seems to be more difficult 
to equilibrate than nearest-neighbor models. For exam- 
ple, we performed 10 7 Monte Carlo sweeps to equilibrate a 
system of around 200 dipoles. To compare, with the case 
of the Heisenberg Edwards- Anderson spin glass, around 
10 7 Monte Carlo sweeps, with both overrelaxation and 
heatbath sweeps counted as a Monte Carlo sweep, were 
used to equilibrate a system of 32,768 spins. 2 '' In simu- 
lations of the Ising Edwards- Anderson spin glass around 
6.5 • 10 6 Monte Carlo sweeps were used to equilibrate a 
system of 8000 spins. 11 Further progress in exploring the 
freezing in Ising 46 ' 53 and Heisenberg (this work) dipolar 
spin glasses will necessitate more sophisticated methods. 
We hope that our present work motivate such develop- 
ments. 
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Appendix A: MAGNETIZATION AND 
STAGGERED MAGNETIZATION 

As system of dipoles placed on the fully occupied SC 
lattice, or when the fraction of vacant sites is sufficiently 
low, orders antiferromagnetically. 100 ' 101 To rule out the 
presence of a long-range order we calculate the magne- 
tization and staggered magnetization. The thermal and 
disorder averaged magnitude of magnetization is defined 
as 



M 



1 N 



(Al) 



where (. . .) denotes thermal averaging and [. . .] is a dis- 
order average. 

The antiferromagnetic ground state (GS) of a system 
of dipoles on fully occupied SC lattice is described by a 
spin vector with the following components: 100 ' 106 



Sf = rf sin 9 cos t 
Sf = rf sin 8 sin < 

Sf = T/COS0, 



(A2) 



Such GS has two global rotational degrees of freedom: 
polar angle, 9, and azimuthal angle, 4>. The sublattice 
and direction indexing vector, r, = [rf , rf , rf], is given 

by 



[(-l) r . ?+r « z , (-ifi+ r t , {-ifi+ r \ 



(A3) 
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Figure 14: (color online). Magnetization, M, and staggered 
magnetization, Af sta g, x=0.0625. 
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Figure 15: (color online). Magnetization, M, and staggered 
magnetization, M sta g, £=0.125. 



where is the position of site i, measured in units of 



lattice constant, and its vector components, 



and 



r\ , on SC lattice, are all integers. The staggered magne- 
tization, which is indicating ordering described by Eqs. 
( A2) , using sublattice and direction indexing vector r, of 
Eq. (A3), is given by 



AL 



stag 



1 N 



(A4) 



In Figs. 14 and 15 we plot the magnetization, M, 
and the staggered magnetization, M s tag, for x=0.0625 
and 2=0.125, respectively. M has a small value that de- 
creases with system size, L. This indicates that nonzero 
magnetization is just a finite-size effect and not an in- 
dication of long-range order. Furthermore, M remains 
constant at all temperatures and does not increase below 
T g . Mstag, similarly to M, decreases with increasing sys- 
tem size, L, and there are no features indicating order- 
ing transition. The fairy large magnetization indicates 
that the finite-size effects are large, and thus the scal- 
ing corrections are expected to be large. The staggered 
magnetization, M sta g, is smaller than the magnetization, 
M. Relatively large magnetization can indicate forma- 
tion of short-range ferromagnetic blocks. Ferromagnetic 
spin blocking has been observed in experimental studies 
of diluted dipolar Heisenberg SG systems Eu^Sri-^S 57 ' 58 
and EuzCai-zBg. 51 ' 



Appendix B: PERIODIC BOUNDARY 
CONDITIONS AND SELF-INTERACTION 



We consider a dipolar Hamiltonian of the form 



n = \e d — V « v sun)^), (Bi) 

where /U and v are vector components, n,v=x, y or z. H 
can be written as 



H = U E ^! s i s h 



(B2) 



or shorter 



(B3) 



where 



£ij — £{ r ij) — 



1,3 



C "ij'iJ 



(B4) 



To impose periodic boundary conditions, we replace the 
interaction matrix, C%j , with 



Ei 8jj \r-jj +n\ - 3(r i:j + n)^(r aj - + n) v 



(B5) 

where n = kLx + ILy + mLz; k, /, m are integers and 
x, y and z are unit vectors. L is the linear dimension 
of the cubic simulation box in units of a, the linear size 
of the cubic unit cell. ^2 n ' means that the summation 



does not include the n=0 term for i = j, where r< 



0. 



One must be aware of the presence of the (n^ 0) self- 
interaction term 



3n»n u 



(B6) 



The self-interaction term describes the interaction of a 
dipole with its own periodic images replicated outside the 
simulation box. For a cubic simulation box it reduces to a 
simple form V?" = LuS^ . To show that the off-diagonal 
terms are zero, for ^J^vwe write 



E 



n^n u + {-n»)n v 



n 



(B7) 



And further, in cubic symmetry, all three directions x, y 
and z are equivalent; hence, Lff = L\f = L'". 
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Appendix C: EWALD SUMMATION 

We wish to calculate the lattice sum 



jliV 



E 



' 6i 



n\ - 3(r y + ny(n 



n 



(CI) 

Again, the prime symbol in the summation sign means 
that for i = j the sum does not include the n=0 term. 
Noting that 



V V I-felZ 

v li v V — r 

r 



3r„rv 



we write 



n 



(C2) 



(C3) 



hence, we may calculate the lattice summation for a 
Coulomb potential and obtain the sums for dipolar in- 
teractions by taking derivatives afterwards. 

The infinite sum (C3) is conditionally convergent, 
meaning that the result depends on the asymptotic or- 
der of summation. The Coulomb or dipolar potential 
is slowly decaying at large distances; hence, with di- 
rect summation, it converges slowly. To alleviate these 
problems, the summation is performed using the method 
introduced by Ewald. 77 ' 78 ' 79 ' 80 In the Ewald technique, 
we separate the summation into two rapidly convergent 
sums: one performed in the direct (real) space and the 
other sum performed in the reciprocal space. Here we 
show only a simplified derivation, rigorous mathematical 
proofs and detailed discussions can be found in Ref. [80]. 

Using the relation 



'dp, 



(C4) 



we write 



_ rV ^ + erfc(ar) 



where 



erfc(a;) = 1 — erf(cc) 



dy 



(C5) 



(C6) 



is the complementary error function. The second term in 
Eq. (C5), for large a, is decreasing fast with increasing 
r; hence, it converges rapidly in the summation over n. 
The first term falls to zero slowly with increasing r, but 
it converges rapidly in a reciprocal space summation for- 
mulation. The splitting parameter a is chosen such that 



both real space and reciprocal space sums are converg- 
ing equivalently rapidly. To obtain the reciprocal space 
summation term we use the relation 



: ^ e -(r+n)V 



27T Y,p~ 3e ~ K2/4p2elK ' r > 



L 3 



(C7) 



K 



where K are the reciprocal lattice vectors, n = L(kx + 
ly + mz); k, I, m are integers and x, y and z are unit vec- 
tors. Some care must be taken to account for the partic- 
ular case of r=0 that corresponds to the self-interaction 
(B6). In that case, n=0 should be excluded from the 
summation (C3) and we write 

_Ly- e _ (r+n) v = -3 e -KV4p 2 ^-r._ 2 e -„v 

Jtt ^— ^ L 3 yhr 



Noting that 



dpp- 3 e~ K2 ^ 2 



K 2 ' 



-K 2 /4a 2 



(C8) 



(C9) 



we can write 



E 



' i 



E 



erfc(a \nj + n\) 



E 

2a . 

--^Oi. 



^ C -K 2 /4a 2 c iK-v, 



L 3 K 



(CIO) 



The divergent, K=0 term in the reciprocal lattice sum- 
mation is omitted. 

To calculate the dipolar sum in (C3), we need to take 
derivative of expression (C8). To start, we compute 



V M V„ 



where 



erfc(ar) _ 8^ v B(r)r 2 - C{r)r ll r v 



B(r) =erfc(r) + ^e- Q2r2 , 

\/7T 



and 



C(r) = 3erfc(r) H = e 



(Cll) 



(C12) 



(C13) 



For the reciprocal space part we compute 

-V^ v e lK ' r = K^K v e iK ' r . 

To obtain the self term (the last term in Eq. C8) we 
write 

- V p X7„e- r2 ? 2 = 2p 2 (6^ - 2p\r v ) e^? 2 , (C14) 
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and integrating (see Eq. C9) we get — -^^S^Sij . Finally, 
we have 



T fiV 



2^ ZE 



47T ^-V K^Ky lf2 tA „2 



L 3 ^ K 2 

4a 3 



'"e-^/^e^ (C16) 



(C17) 



Similarly to the Coulomb case (CIO), the divergent -KT=0 
term is omitted in the reciprocal space summation. 

The effect of the magnetic polarization of the surface 
does not vanish in the thermodynamic limit. To model 
the experimental case of a spherical sample, a direct (real 
space) sum (CI) can be computed via summing over se- 
ries of spherical shells of radius r k , where each shell con- 
sist of all vectors n such that r k < \n\ < r k +i. In the 
Ewald method, to obtain a result equivalent to such a 
summation, the surface contribution to the total energy 
should be included, and it is of the form 79 



[/(surf) 



2?r 



(2e' + l)L- 



E 



where e' is the magnetic permeability of the surround- 
ing medium. In the case of a long cylindrical shape 
the surface term is zero. In our simulations we set the 
surface term to zero and are therefore implicitly con- 
sidering a long cylindrical sample. In practice, we set 
e' = co, infinite magnetic permeability outside the con- 
sidered system, the so-called "metallic boundary condi- 
tions", by analogy to the physical situation with electric, 
as opposed to magnetic dipoles. 



Appendix D: OVERRELAXATION 

It has been reported that supplementing canonical 
Metropolis spin updates with computationally inexpen- 
sive "overrelaxation" steps of zero energy change can 
substantially reduce autocorrelation times. M,,>il Unfortu- 
nately, this technique does not provide much of a perfor- 
mance improvement in the case of long-range interactions 
and cannot be used when periodic boundary condition 
are imposed on a system characterized by dipolar inter- 
action with non-cubic lattice symmetry. 

In overrelaxation update, a new spin direction, S^, is 
obtained by performing a reflection of the spin at site i, 
Si, around the local dipolar field vector, Hi, 



S> = -S i + 2?^H i . 
The local dipolar field is given by Eq. (5), 



H k - ^2L kj Sj, 

3^k 



(Dl) 



(D2) 



where the tensor L k j stands for the dipolar interaction, 
as defined in Eq. (2). Using Eq. (D2), the finite-size 
Hamiltonian of Eq. (3) can be written in the form 



H = -- ^2 Sk ' Hk ~ o E Sk ^ kk ' 



(D3) 



where the local dipolar field, H k , does not include the 
self-term and the self term is written explicitly. Let us 
consider an overrelaxation move of spin Si. To make the 
effect of the spin move clear, we write the energy of a 
spin configuration before the spin move in the form 



E 



-Si ■ H t - \ J2 Sk ■ Hk ~ \ E SklkkS ®±) 



k^i 



which is just Eq. (D3) rewritten with the term for spin 
Si excluded from the summation and written explicitly. 
After changing spin Si to S 1 ,-, according to Eq. (Dl), we 
have 



E' = —SI ■ ^ - S k ■ H' k - \ S' k L kk S' k , (D5) 

k^i k 

where H' k are updated dipolar fields; H[ = Hi and for 

H'^Hk + LkiiSi-Si). (D6) 
Combined together, Eqs. (Dl), (D4), (D5) and (D6) give 



E' — E — — yS'iLiiSl — SiLuSi 



The energy does not change only if .S'/.,,.S'' Sil-aS 
This is the case when, for each [i,v, = 



(D7) 

i = 0. 
or for 



diagonal La, = LaS^, (recalling = \Si\ = 1), 
which is satisfied in the case of cubic lattice symmetry 
(see Appendix B). 

The fact that we do not use the overrelaxation method 
does not cause a large decrease of efficiency in our simula- 
tion. In the case of long-range interaction, the reflection 
(Dl) would have to be followed by the recalculation of 
dipolar field, H' k , of Eq. (D6). A similar lattice sum has 
to be performed in the case of Metropolis updates. Most 
of the computation time is spent on doing such lattice 
sums. Hence, even if it was doable, an overrelaxation 
move would be practically as computationally expensive 
as a Metropolis update. 



Appendix E: HEATBATH ALGORITHM 

In the original Metropolis algorithm a random config- 
uration update is attempted and it is accepted with a 
probability that depends on the change of energy follow- 
ing such configuration change. The updates lowering the 
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energy are always accepted while, if the energy is to in- 
crease, the acceptance probability is 



P(AE) = exp(-/3A£0, 



(El) 



where (3 = 1/k^T. The probability exponentially de- 
creases with an increase of the energy change, AE. Thus, 
to obtain a sufficient acceptance rate, the attempted 
moves have to be sufficiently small. Usually the con- 
figuration update is chosen in such way that the accep- 
tance rate is close to 50%. In many applications a better 
way of performing a local spin updates is the "heatbath" 
algorithm, 88 ' 89 ' 90 where the new direction of a spin is 
drawn from a suitable probability distribution such that 
the new configuration energy is distributed according to a 
Boltzmann weight. In the case of isotropic (0(3)) Heisen- 
berg model, the distribution of angle 9 between the local 
dipolar field, Hi, and the spin vector Si can be calcu- 
lated analytically. 88 ' 89 ' 90 In such an isotropic case, the 
Hamiltonian can be written as 



H 



(E2) 



where Hi is the interaction field and Hi does not depend 
on spin Si. We did not include here the self- interaction 
term because the calculation shown below is not possible 
with a general self-interaction term included. The case 
of diagonal self-interaction term, as in the case of cubic 
symmetry, will be discussed at the end of this appendix. 
In the case of long-range interaction we write 



Hi — Lij Sj 



(E3) 



It is convenient to describe spin S, in polar coordinates, 
9 and <j>, with the polar axis along the local dipolar field, 
Hi. The energy of spin i in the field of other spins is 



Ei = S, H, = -flicos(fl), 



(E4) 



where 9 is the polar angle defined as the angle between 
Si and Hi. We wish to randomly choose Si such that 
the probability distribution of the energy (E4) given by 
Boltzmann distribution. The energy does not depend on 
the azimuthal angle, (j>; hence, <j) is randomly chosen from 
the uniform distribution on the interval [0, 2ir]. The polar 
angle, 9, is chosen such that x = cos{0) is given by the 
probability distribution 



P(x) 



f^dxeP 11 ^ 2smh/3iJ i 



(E5) 



where j3 = 1/T. To obtain random variable x drawn from 
distribution (E5), we calculate the cumulative distribu- 
tion 



F{x) 



P(x')dx' 



e l3HiX _ e -f3H> 



e f3H* - e - 

and we reverse r = F(x), where r <E [0, 1] is a uniformly 
distributed random number. We obtain 88 ' 89 ' 90 



-0Hi 



(E6) 



1 



In [1 + r (e 



2f3Hi 



1 



1. 



(E7) 



Having chosen 4> and 9, we need to compute the vector 
components of the spin and rotate them to the global co- 
ordinate system. Using the coordinates <f> and 9, relative 
to the local molecular field, H^, we compute the new 
spin vector, Si = (S x , S y , S z ), as follows. Let 4>h and 9h 
denote azimuthal and polar angle of vector Hi in global 
coordinates. A possible choice of the local coordinates 
x! , y' and z' , having i' axis along Hi is 



x' = 

y = 



z' 



cos(9h) cos(4>h)x + cos(9h) sin(^)y — sin(6*#).z 
= - sin((f> H )x + cos((/)H)y (E8) 
= s\u(9h) cos(4>h)x + sin(##) sin(^)y + cos(9h)z. 



The new spin, Si, in local coordinates is 

Si = sm{9) cos(4>)x' + sin(0) sin(0)y' + cos(9)z', (E9) 

and finally, combining Eq. (E8) and Eq. (E9), we get 

S x = 9 C08((f)ff) — shi(9) sin(0) sin(^>#-), 

S y = 8sin(^jj) +sin(0)sin(0)cos(^ J? ), (E10) 

S z = — sm(9) cos(4>) sin(^H) + cos(9) cos(6h), 

where 

9 = sin(0) cos(0) cos(9 H ) + cos(fl) sin(6» H ). (Ell) 

Hamiltonian (E2) with dipolar field (E3) does not in- 
clude a self-interaction term. In the case of dipolar in- 
teraction with periodic boundary condition we have to 
include such a self-interaction term as in the Hamilto- 
nian of Eq. (D3). For clarity, we write this Hamiltonian 
again 



H = - 



2 ^ 



Si Hi + S, l.,,S, 



(E12) 



For Hamiltonian (E12), as opposed to Hamiltonian (E2), 
for a general form of the matrix Lu , the cumulative dis- 
tribution (E6) cannot be integrated and reversed analyt- 
ically. However, for cubic symmetry, the self-interaction 
term is of the form L^" 
to 



Lu5^ v and Eq. (E12) reduces 



n 



Li 



(E13) 



hence, it is of the form (E2), with just a constant inde- 
pendent of spin configuration added, and the heatbath 
method can be applied. 

Although we have cubic symmetry in our system and 
the heatbath algorithm could in principle be used, we 
decided to employ a method generally applicable for 
any lattice symmetry and we did not use the heatbath 
method. 
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